Image analysis

ABSTRACT

A method for the automated analysis of digital images, particularly for the purpose of assessing mitotic activity from images of histological slides for prognostication of breast cancer. The method includes the steps of identifying the locations of objects within the image which have intensity and size characteristics consistent with mitotic epithelial cell nuclei, taking the darkest 10% of those objects, deriving contours indicating their boundary shape, and smoothing and measuring the curvature around the boundaries using a Probability Density Association Filter (PDAF). The PDAF output is used to compute a measure of any concavity of the boundary—a good indicator of mitosis. Objects are finally classified as representing mitotic nuclei or not, as a function of boundary concavity and mean intensity, by use of a Fisher classifier trained on known examples. Other uses for the method could include the analysis of images of soil samples containing certain types of seeds or other particles.

FIELD OF THE INVENTION

The present invention relates to the automated analysis of digitalimages. It is more particularly concerned with the automatedidentification of mitotic activity in digital images of histological orcytology specimens and most particularly for the purpose of assessingthe presence and severity of cancer in breast tissue, and it is in thiscontext that the invention is principally described herein. Theinvention may, however, also find application in the assessment of otherforms of cancer, such as colon and cervical cancer, and in the analysisof various other kinds of structure presenting image components whichare amenable to identification in a similar way, for example in theanalysis of soil samples containing certain types of seeds or otherparticles.

BACKGROUND AND SUMMARY OF THE INVENTION

Many thousands of women die needlessly each year from breast cancer, acancer from which there is theoretically a high probability of survivalif detected sufficiently early. If the presence of cancerous tissue ismissed in a sample, then, by the time the next test is undertaken, thecancer may have progressed and the chance of survival significantlyreduced. The importance of detecting cancerous tissue in the samples cantherefore not be over-emphasised.

A typical national breast screening programme uses mammography for theearly detection of impalpable lesions. Once a lesion indicative ofbreast cancer is detected, then tissue samples are taken and examined bya trained histopathologist to establish a diagnosis and prognosis. Moreparticularly, one of the principal prognostic factors for breast canceris the extent of mitotic activity, that is to say the degree ofepithelial cell division that is taking place. A histopathological slideis effectively a “snapshot” representing a very short time interval in acell division process, so the chance of a particular slide showing aparticular phase of mitotic activity is very small; if such a phase isin fact present in a slide, that is a good indicator of how fast apotential tumour is growing.

In the existing manual procedure for scoring mitotic activity ahistopathologist places a slide under a microscope and examines a regionof it (referred to as a tile) at a magnification of ×40 for indicationsof mitoses. Typically ten different tiles from the tissue sample areexamined and a total count is made of the number of cell divisionswhich, in the histopathologist's opinion, are seen to be taking place inthe ten tiles. This is then converted to an indication of cancer gradetypically in accordance with the following table: Number of MitoticCells per Ten Tiles Cancer Grade 0 to N Grade 1 (N + 1) to M Grade 2 >MGrade 3where Grade 1 is the least serious and Grade 3 is the most serious.Values of N and M are typically 5 and 10 but will vary in differentschemes depending on the size of the tiles being observed.

This is, however, a time consuming, labour intensive and expensiveprocess. Qualification to perform such examination is not easy to obtainand requires frequent review. The examination itself requires theinterpretation of colour images by eye, a highly subjective processcharacterised by considerable variations in both inter, andintra-observer analysis, i.e. variances in observation may occur for thesame sample by different histopathologists, and by the samehistopathologist at different times. For example, studies have shownthat two different histopathologists examining the same ten samples maygive different opinions on three of them, an error of 30%. This problemis exacerbated by the complexity of some samples, especially in marginalcases where there may not be a definitive conclusion. If sufficienttrained staff are not available this impacts upon pressures to completethe analysis, potentially leading to erroneous assessments and delays indiagnosis.

These problems mean that there are practical limitations on the extentand effectiveness of screening for breast cancer with the consequencethat some women are not being correctly identified as having the diseaseand, on some occasions, this failure may result in premature death.Conversely, others are being incorrectly diagnosed with breast cancerand are therefore undergoing potentially traumatic treatmentunnecessarily.

It is thus an aim of the invention to provide an automated method ofimage analysis which can be embodied in a robust, objective andcost-effective tool to assist in the diagnosis and prognosis of breastcancer, although as previously indicated the invention may also findapplication in other fields.

In one aspect the invention accordingly resides in a method for theautomated analysis of a digital image comprising an array of pixels,including the steps of: identifying the locations of objects within theimage which have specified intensity and size characteristics; definingregions of specified extent within the image which contain respectivesaid objects; deriving from the data within respective said regions oneor more respective closed contours comprising points of equalintensities; and estimating the curvature of at least one respectivesaid contour within respective said regions at least to produce ameasure of any concavity thereof.

As will be understood from the ensuing detailed description of apreferred embodiment, such a method is of use in identifying mitoticcell nuclei in digital images of histopathological slides.

The invention also resides in apparatus for the automated analysis of adigital image comprising means to perform the foregoing method and in acomputer program product comprising a computer readable medium havingthereon computer program code means adapted to cause a computer toexecute the foregoing method and in a computer program comprisinginstructions so to do.

These and other aspects of the invention will now be more particularlydescribed, by way of example, with reference to the accompanyingdrawings and in the context of an automated system for grading cancer onthe basis of the numbers of mitotic epithelial cell nuclei in digitalimages of histopathological slides of potential carcinomas of thebreast.

BRIEF DESCRIPTION OF THE DRAWINGS

In the drawings:

FIG. 1 is a block diagram of an automated process in accordance with theinvention for measuring mitotic activity for patient diagnosis;

FIG. 2 is a more detailed block diagram of the main stages in themitosis detection and measurement block of FIG. 1;

FIGS. 3 and 4 are simplified visualisations of the contour selectionstage of the process of FIG. 2; and

FIG. 5 illustrates a decision boundary in a Fisher classifier as used ina later stage of the process of FIG. 2.

DETAILED DESCRIPTION

FIG. 1 shows a process for the assessment of tissue samples in the formof histopathological slides of potential carcinomas of the breast. Theprocess measures mitotic activity of epithelial cells to produce aparameter for use by a pathologist as the basis for assessing patientdiagnosis. It employs a database 1, which maintains digitised image dataobtained from histological slides. Sections are cut from breast tissuesamples (biopsies), placed on respective slides and stained using thestaining agent Haematoxylin & Eosin (H&E), which is a common stain fordelineating tissue and cellular structure.

To obtain the digitised image data for analysis, a histopathologistscans a slide under a microscope and at 40× magnification selectsregions of the slide which appear to be most promising in terms ofanalysing mitotic activity. Each of these regions is then photographedusing the microscope and a digital camera. In one example a ZeissAxioskop microscope has been used with a Jenoptiks Progres 3012 digitalcamera. This produces for each region a respective digitised image inthree colours, i.e. red, green and blue (R, G & B). Respective intensityvalues in the R, G and B image planes are thus obtained for each pixelin an array. In the preferred embodiment there is an eight-bit range ofpixel intensities (values of 0 to 255) for each colour and the arraycomprises 1476 pixels across by 1160 pixels down, with a pixel size of220 nm square. The image data is stored temporarily at 1 for later use.Ten digitised images (electronic equivalents of tiles) are required forthe detection and measurement of mitotic activity at 2, which thenprovides input to a diagnostic report at 3. In principle the processingstages to be described in detail below can operate on a single waveband(R, G or B) image or a combination of them. In practice, however, thered waveband has been found to contain the most information fordiscriminating between mitotic and other cells when stained with H&E,and is assumed to be used in the following description. The process canbe performed in a suitably programmed personal computer (PC) or othergeneral purpose computer of suitable processing power or in dedicatedhardware.

FIG. 2 shows in more detail the processing stages comprised in the block2 of FIG. 1. They are carried out for each of the ten digitised imagesreferred to above and will now be described for one such image (file).To aid in the understanding of this description it is recalled that theaim of the process is to identify and count the number of mitoticepithelial cell nuclei (if any) in each tile. In images acquired asdescribed above such nuclei generally appear darker than normalepithelial cell nuclei, and also have a different shape. Normal nucleiare generally convex with smooth boundaries while mitotic nuclei aremore irregular in shape and have ragged boundaries. However, it is notalways the case that mitotic epithelial cell nuclei are the darkestobjects in a given file; for example stromal cells, lymphocytes andnecrotic cells may be darker. Given the relatively low numbers ofmitoses which may be present in any given tile and yet may indicateserious disease it is important that as many as possible are correctlyidentified while at the same time minimising the number of any normalcell nuclei or other objects incorrectly identified as mitotic.

Location of Candidate Cell Nuclei

Referring now to FIG. 2, the first processing stage 21 consists oflocating all possible candidate cell nuclei. The approach adopted foridentifying the locations of potential mitotic nuclei is based on thefact that they are generally darker than average nuclei. Mitotic nucleiappear in the image as solid dark objects (i.e. dark all the waythrough) most of the time, or instead occasionally they form groups ofsmall dark clumps. Hence the aim is to find concentrations of darkpixels; these are not necessarily connected groups of dark pixels, but aregion containing a sufficient number of clustered dark pixels.

There are various methods of doing this. One example is simplegrey-level segmentation, where a threshold is chosen and only thosepixels having grey-levels below this threshold are selected. Thedrawback of this approach is that some mitotic nuclei are notparticularly dark, but are only distinguishable from their shapecharacteristics. Choosing a threshold sufficiently low to detect suchnuclei would yield an excess of clutter.

The preferred approach is to use the multiresolution blob filteringdescribed below. However, as will be apparent to those skilled in theimage processing art, the present invention may be practised withoutemploying this particular technique. Alternatives include the processesdescribed as mitotic cueing in our copending United Kingdom patentapplication no. 0226787.0. The general principle is, given that theapproximate size of the nuclei is known, to apply a radially-symmetricfilter whose output is large in magnitude when there is a region oflocal brightness or darkness whose shape and size approximately matchesthat of the filter. This filter should be a difference filter with zeromean, so areas of constant intensity are suppressed.

The method will now be described in terms of a specific implementation,namely a multi-scale blob filter as known e.g. from “Multiresolutionanalysis of remotely sensed imagery”, J. G. Jones, R. W. Thomas, P. G.Earwicker, Int J. Remote Sensing, 1991, Vol 12, No 1, pp 107-124. Theprocess will be described for filtering the image using a particularsize of blob filter, where these are defined at successive octave(powers of 2) scale sizes.

The recursive construction process for the multi-scale blob filterinvolves two filters; a 3×3 Laplacian blob filter (L) and a 3×3smoothing filter (s) as defined below. ${L = \begin{pmatrix}{- 1} & {- 1} & {- 1} \\{- 1} & 8 & {- 1} \\{- 1} & {- 1} & {- 1}\end{pmatrix}},{s = {\frac{1}{16}\begin{pmatrix}1 & 2 & 1 \\2 & 4 & 2 \\1 & 2 & 1\end{pmatrix}}}$

These two filters form a basis for filtering over a set of octave scalesizes, according to the following process:

To enhance blob-shaped objects at the original resolution (octave scale1, pixel size 1), the image is correlated with the 3×3 Laplacian filteralone:${F_{1}\left( {m,n} \right)} = {\sum\limits_{i = {- 1}}^{I}{\sum\limits_{j = {- 1}}^{I}{{I\left( {{m + i},{n + j}} \right)}*{L\left( {{i + 2},{j + 2}} \right)}}}}$where I is the original image, and the range of the indices m and n isset such that the indices in the summation above are always within theoriginal image dimensions (so m and n start at 2). Values of thefiltered image at locations outside these ranges are set to zero.

For computational efficiency, multiplications by ±1 need not beperformed explicitly. Thus the filter output value for a pixel locatedat position (i,j) is given by:F ₁(i,j)=8.I(i,j)−I(i−1,j−1)−I(i−1,j)−I(i−1, j+1)−I(j−1)−I(i,j+1)−I(i301,j−1)−I(i+1,f)−I(i+1,j+1)

To enhance blob-shaped objects at a resolution one octave above theoriginal (octave scale 2, pixel size 3), the image is first correlatedwith the 3×3 smoothing filter (s), forming a smoothed image (S₂). The3×3 Laplacian blob filter (L) is then expanded by a factor of two, bypadding it with zeros, to form a 5×5 filter [−1 0 −1 0 −1; 0 0 0 0 0; −10 8 0 −1; 0 0 0 0 0; −1 0 −1 0 −1]. This is then correlated with thesmoothed image (S₂) to form a filtered image (F₂), but for computationalefficiency, only the non-zero filter coefficients are used, thus:${S_{2}\left( {m,n} \right)} = {\sum\limits_{i = {- 1}}^{I}{\sum\limits_{j = {- 1}}^{I}{{I\left( {{m + i},{n + j}} \right)}*{s\left( {{i + 2},{j + 2}} \right)}}}}$${F_{2}\left( {m,n} \right)} = {\sum\limits_{i = {- 1}}^{I}{\sum\limits_{j = {- 1}}^{I}{{S_{2}\left( {{m + {2i}},{n + {2j}}} \right)}*{L\left( {{i + 2},{j + 2}} \right)}}}}$where I is the original image, and the range of the indices m and n isset such that the indices in the summation above are always within theoriginal image dimensions (so m and n start at 4). Values of thefiltered image at locations outside these ranges are set to zero.

The above double correlation is equivalent to a single correlation ofthe original image with a 7×7 filter formed from correlating theexpanded 5×5 Laplacian with the 3×3 smoothing filter, but this largerfilter is never formed explicitly.

To enhance blob-shaped objects at a resolution two octaves above theoriginal (scale 3, pixel size 7), the smoothing filter is expanded by afactor of 2 in the same manner as the Laplacian above, then correlatedwith the smoothed image (S₂) above to give a lower-resolution smoothedimage (S₃) thus:${S_{3}\left( {m,n} \right)} = {\sum\limits_{i = {- 1}}^{I}{\sum\limits_{j = {- 1}}^{I}{{S_{2}\left( {{m + {2i}},{n + {2j}}} \right)}*{s\left( {{i + 2},{j + 2}} \right)}}}}$

Following this, the 5×5 Laplacian filter is expanded by a factor of 2 bypadding with zeros to form a 9×9 filter, which is correlated with thesmoothed image (S₃) in the same computationally efficient manner, thus:${F_{3}\left( {m,n} \right)} = {\sum\limits_{i = {- 1}}^{I}{\sum\limits_{j = {- 1}}^{I}{{S_{3}\left( {{m + {4i}},{n + {4j}}} \right)}*{L\left( {{i + 2},{j + 2}} \right)}}}}$

This process is repeated to obtain results at successive octave scales,namely expanding .both the smoothing filter and the Laplacian blobfilter each time.

The above process may be used to produce a “blob-filtered” image at anyof the required octave scales. Objects of interest (i.e. clusters ofdark pixels in this case) will have the greatest values of the magnitudeof the filter output The locally strongest filter output will occur atthe centre of the object of interest. Individual objects (called“blobs”) are now identified by finding local minima of the filteroutput, where each blob is assigned a position and intensity, the latterbeing the value of the filter output at the local minimum.

In this application, objects which are dark relative to the backgroundare identified by finding local minima of the filter output at onechosen scale, in this instance octave scale 5 (a nuclear size of 31pixels across).

For computational efficiency, the spatial resolution of the image isreduced prior to blob filtering, using the following thinning method.Each reduction in resolution by a factor of two (a “thin”) is achievedby firstly correlating the image with the 3×3 smoothing filter(s), thensub-sampling by a factor of two. The formula for a single thin is:${T\left( {m,n} \right)} = {\sum\limits_{i = {- 1}}^{I}{\sum\limits_{j = {- 1}}^{I}{{I\left( {{{2m} + i},{{2n} + j}} \right)}*{s\left( {{i + 2},{j + 2}} \right)}}}}$where I is the original image, T is the thinned image, and the indices mand n range from 1 to the dimensions of the thinned image. Eachdimension of the thinned image is given by subtracting 1 or 2 from thecorresponding dimension of the original image, depending on whether thisis odd or even respectively, then dividing by 2.

In this instance, the image is reduced in resolution by a factor of 4,by applying the above process twice, firstly on the original image, andthen again on the resulting thinned image, to produce a new image whoselinear dimensions are a quarter of the size of the original (area is1/16). For example, an original image of a tile of size 1476×1160 pixelswould become 368×289 pixels. Blobs are now extracted as described abovefrom the reduced-resolution image at octave scale 3 (7×7 pixels across),this being equivalent to extracting scale 5 blobs from the originalimage, but being more computationally efficient.

This process identifies all objects which form dark clusters of therequisite size, which may include not only mitotic epithelial cellnuclei but also background clutter, stromal cells, lymphocytes and/ornecrotic cells, plus fainter normal epithelial nuclei which are not ofinterest. Since it is known that the mitotic nuclei are likely to bemuch darker than average, only the darkest 10% of blobs are selected forfurther analysis. This is achieved by sorting the blobs into ascendingorder of filter output (so the darkest occur first), using a QuickSortalgorithm (such as described in Klette R., Zamperoniu P., ‘Handbook ofImage Processing Operators’, John Wiley & Sons, 1996), finding the10^(th) percentile of the sorted values, and choosing all blobs darkerthan this percentile.

Segmentation and First Clutter Rejection

The next processing stage 22 aims to find an approximate segmentation ofthe image, to separate regions (defined as connected sets of pixels)potentially associated with the cell nuclei of interest, from thebackground. The full-sized original image (red component) is used at thecommencement of this stage.

Firstly, a grey-level threshold is selected. This is achieved bychoosing a set of 15×15 pixel neighbourhoods centred on each of theblobs selected at the end of stage 21, collating all pixels within allthese neighbourhoods into a single list, and computing the meangrey-level of these pixels.

A new thresholded binary image is now produced. Pixels in the redcomponent of the original image whose grey levels are below (darkerthan) the threshold mean computed above are set to 1; remaining pixelsare set to 0.

Connected component labelling is now applied to this binary image. Thisis a known image processing technique (such as described in A Rosenfeldand A C Kak, ‘Digital Picture Processing’, Vols. 1 & 2, Academic Press,New York, 1982) which gives numerical labels to connected regions in thebinary image, these being groups of connected pixels whose values areall 1. An 8-connectedness rule is used, so pixels are deemed to beconnected when they are horizontally, vertically, or diagonallyadjacent. Each region corresponding to a selected blob from stage 21 isassigned a separate label, enabling pixels belonging to those regions tobe identified. The following region properties are then computed:

-   -   Area= number of pixels within the region    -   Thickness=minimum thickness of the region, defined thus: for        each pixel in the region, find the minimum distance from that        pixel to the outside of the region. Thickness is then defined to        be the maximum of these minimum distances. Note that thickness        is not the same as width; for a rectangle the thickness is half        the width, and for a circle the thickness is the radius.

Regions whose area is less than 190 pixels or whose thickness is lessthan 4 pixels are rejected, these being too small to be mitotic cellnuclei.

At this stage the mean grey-level of the pixels within each region isalso calculated from the red component of the original image. Theoverall mean and standard deviation of these mean grey-levels is thenfound for later use in grey-level normalisation (stage 25).

Contour Selection

The next processing stage 23 incorporates two levels of contourselection to gain a better representation of the actual shape of theboundary of each remaining object at both low and high resolutions.Firstly, a low-resolution (large-scale) contour is computed, which givesan approximate shape, and secondly a high-resolution (small-scale)contour is found which gives a more accurate boundary representation.Following consistency checks between the two contours, attributes of theboundary are then measured from the small-scale contour.

For each of the objects remaining after stage 22, a local region ofinterest (ROI) is selected. This ROI is centred on the nominal centre ofthe object (as found in stage 21), and has an extent of 50 pixels ineach direction, the region size being truncated as necessary to ensurethe ROI lies within the bounds of the original image. This allows ROlswhich would otherwise overlap the edges of the image to be included.Alternatively the ROls could be defined by taking the regions identifiedin stage 22 and adding a border of a selected number of pixels. Ineither case, it is desirable that the ROls exceed the size of thoseregions somewhat in order to ensure the generation of the low-resolutioncontours.

To find a low-resolution representation for the boundary of each object,the region defined by the ROI above is used to define a sub-image withinthe output of the blob filter (stage 21). This sub-image will consist ofboth positive and negative grey levels. Contours at two levels withinthis sub-image are then sought, namely at levels 0 and −10 which havebeen found to be best experimentally. By virtue of the operation of theblob filter in stage 21, the zero level contour in the respectivesub-image is that contour which exhibits the highest edge strength. Acontour is a curve consisting of points of equal value for some givenfunction; in this case the function is defined by the grey-level pixelvalues. In this embodiment, the Matlab® contour function is employed butany contouring algorithm can be used which returns contours in the sameform, as a set of contiguous points ordered around the contour; (Matlab®is a well known computational tool from The MathWorks, Inc.). Matlab®returns a set of locations with sub-pixel resolution which are in orderof location around the contour, i.e. traversing the set of locations isequivalent to walking around the contour. Contours are only treated asvalid if they satisfy all the following four conditions:

-   -   they form closed loops within the ROI, i.e. the last contour        point is the same as the first contour point;    -   they are consistent with the location of the object (there is at        least one contour point whose distance from the nominal centre        of the object is less than or equal to 30 pixels);    -   they have a sufficiently similar area to the “nominal area”        found from the grey-level segmentation computed in stage 22 (the        definition of the area within a contour is given later in this        section). The contour area must be at least 50% of the nominal        area;    -   they have the correct grey-level orientation, namely pixels        within the contour are darker than those outside the contour.

The object is retained for further analysis (maintained in list incomputer) only if a valid contour is found from at least one of the twocontour levels (0 and −10). If both contour levels yield a validcontour, then the latter one (−10) is chosen for further use.

To find a high-resolution representation, the region defined by the ROIabove is taken out from the red component of the original image to forma sub-image. Contours are not extracted from the image at its originalresolution, because these have been found to be too rough. Instead, theresulting sub-image is expanded in size by a factor of two, using theMatlab® bilinear interpolation function, to give additional resolution.In bilinear interpolation, to find the values of a selected point not onthe original image grid, its four nearest grid points are located, andthe relative distances from the selected point to each of its fourneighbours are computed. These distances are used to provide a weightedaverage for the grey-level value at the selected point.

This interpolated image is then smoothed before contouring, bycorrelating (see earlier description with reference to stage 21) with a3×3 smoothing filter(s) defined thus: $s = {\frac{1}{16}\begin{pmatrix}1 & 2 & 1 \\2 & 4 & 2 \\1 & 2 & 1\end{pmatrix}}$

Valid contours are then sought at each of several threshold levels. Therange of threshold levels starts at the minimum grey-level within thesub-image, and increases up to the maximum grey-level in steps of 10grey levels, so the actual contour levels are set adaptively. Validcontours are defined in the same manner as for the low-resolutionboundary above.

Having found a set of valid contours at each threshold level, the edgestrength at each point on the contour is estimated. The edge strength ateach image pixel is defined as the modulus of the vector gradient of theoriginal red component of the image I, where the vector gradient isdefined asgra${I = \left( {\frac{\partial I}{\partial x},\frac{\partial I}{\partial y}} \right)},$where the two partial derivatives are obtained from taking differencesin pixel grey-level values in the X and Y directions respectively. Theedge strength at contour points which lie between image pixels isestimated using bilinear interpolation (as described above) from thenearest pixel values. The mean edge strength around the contour is thencomputed. The contour having the greatest edge strength is chosen asbeing the most representative of the boundary of the object. If no validcontours are found, the object is rejected.

A consistency check between the low-resoluton and high-resolutioncontours is then performed. The area within each contour is computedfrom the boundary contour using Green's Theorem (such as described inChap 6 in Vector Analysis and Cartesian Tensors', D. E. Boume and P. C.Kendal, 2^(nd) ed, Nelson, 1977). This gives the following formula forarea:$A = {{{x_{n}\left( {y_{1} - y_{n}} \right)} + {\sum\limits_{i = 1}^{n - 1}{x_{i}\left( {y_{i + 1} - y_{i}} \right)}}}}$where (x_(i),y_(i)) are the contour points.

These areas are then subjected to the following tests:

-   -   High-resolution area>0.6*low-resolution area    -   High-resolution area<1.4*low-resolution area

If either of these tests fail, the object is rejected.

Finally, there is a check that the low- and high-resolution contours foreach object overlap sufficiently. For each object, two binary images areformed. The first binary image is formed by setting the value of pixelswhich lie within the low-resolution. contour to 1, and those pixelsoutside the contour to 0, using the Matlab® function roipoly. The secondbinary image is formed in the same way from the high-resolution contour.The absolute difference of these two images is taken, resulting inanother binary image in which pixels are set to 1 if and only if theylie within one of the contours and not the other, and 0 otherwise.Connected component labelling (see stage 22 description) is applied tothis new binary image to identify separate regions. The thickness ofeach of these regions is computed (as in stage 22). If any regionthickness exceeds 5 pixels, the corresponding object is rejected.

Simplified visualisations of the effects of the stage 23 processing areshown in FIGS. 3 and 4.

In FIG. 3(a) a local region of interest 30 is defined around the nominalcentre 31 of an object 32 which is represented in this Figure by aseries of contours. A second object 33 also appears in the samesub-image. FIG. 3(b) illustrates the low-resolution boundary contour 34computed for the object 32 and FIG. 3(c) the high-resolution boundarycontour 35 computed for the same. FIG. 3(d) illustrates the overlapbetween these two resolutions with the region of difference 36 shaded.In this case the areas within the contours 34 and 35 are sufficientlysimilar and the thickness of the region 36 is sufficiently small to passall the above-mentioned consistency checks and the object 32 will beretained. In effect these checks are showing that the object issufficiently uniformly dark to potentially represent a mitotic cellnucleus. The object 33 will not be treated as valid for the ROI 30because its contours are not consistent with the centre 31. It will,however, be separately analysed within a separate ROI (not shown)defined around its own nominal centre.

In FIG. 4(a) there is another example of a local region of interest 40defined around an object 41. FIG. 4(b) illustrates the (ow-resolutionboundary contour 42 computed for this object and FIG. 4(c) thehigh-resolution boundary contour 43. In this case the area of thecontour 43 is substantially less than that of the contour 42(approximately 0.4) so it fails the first of the above-mentionedconsistency checks and the object 41 will not be processed further. Thisis indicative of a normal epithelial cell nucleus which has a relativelydarker nucleolus (chosen as the high-resolution boundary because of itshigh edge strength) surrounded by a less dark region (chosen as thelow-resolution boundary).

Boundary Tracking

The next processing stage 24 applies a tracking algorithm to thehigh-resolution contour representing the object's boundary for eachobject retained from the previous stage 23 in order to estimatecurvature. The aim is to smooth the boundary and then measure curvature,because simply calculating curvature from the contour segments gives toorough a measurement. For the identification of mitotic cell nuclei thedegree of non-convexity of the boundary is of interest, so the laftermethod of calculation is inappropriate.

The particular algorithm which has been used in the preferred embodimentis based on a Probability Density Association Filter (PDAF), such asdescribed in Y.Bar-Shalom and T. E. Fortmann, “Tracking and DataAssociation”, Mathematics in Science and Engineering series, vol 179,Orlando, Fla., Academic Press, 1988. This type of tracking algorithm isdesigned to estimate the parameters of a chosen object (target) in thepresence of other measurements which are not related to the object inquestion (noise and clutter). In this case the ‘target’ state variablesare the position, orientation and curvature of the object's boundary,and the measurements are the positions of the contour points and theorientations of the lines joining each pair of contour points.

The PDAF filter requires a model for the dynamics of the boundary state.The boundary dynamics is given by a constant curvature (the radius ofcurvature is set to 10 pixels) plus an assumed random perturbation knownas system noise. This noise is determined by the variance of curvatureperturbation, which is chosen according to how irregular the boundary ofa mitotic cell nucleus is expected to be. In the preferred embodimentthe curvature variance is 9 for position and 0.09 for angle (inradians).

As a starting point, it is assumed that for each object potentiallyrepresenting a mitotic cell nucleus a connected set of edge features hasbeen extracted from the image. In this case, edge features are linesegments joining two adjacent contour points. Each edge feature has thefollowing measurements that were made as part of the contour extractionprocess:

-   -   Position x_(m), y_(m) (horizontal and vertical image        coordinates) of the centre of the edge    -   Orientation θ_(m), i.e. the angle between the edge and the        horizontal.

The purpose of the tracker is to estimate the most likely true location,orientation and curvature x_(s), y_(s), θ_(s),κ of the boundary at eachpoint from the above measurements, given that there are measurementerrors with an assumed Gaussian distribution. The following informationvectors are defined:

-   -   The measurement vector z=(x_(m), y_(m), θ_(m));    -   The system state vector x=(x_(s), y_(s), θ_(s),κ).

To use the PDAF filter to do this, the following information about thetrue boundary and the measurement process is required:

-   -   The relationship between the position, orientation and        curvature-of neighbouring points on the boundary (the system        model). This incorporates a transition matrix Φ linking        neighbouring states x and a system noise model that adds extra        random perturbations to x.    -   The relationship between the measurement vector z and the system        state x. This incorporates a transition matrix H linking x to z        and a measurement noise model that adds extra random        perturbations to z.    -   It is assumed that not all of the edge features are associated        with the nuclear boundary; the ones that are not are denoted        clutter.

In its most general form the PDAF processes several measurements z ateach step in estimating x. In this case only one edge feature isprocessed at a time, so there are only two hypotheses to be tested;either the feature is from clutter or from the real nuclear boundary.

The system transition matrix Φ is based on constant curvature, so topredict a neighbouring system state the unique circle or straight linewith curvature κ, tangent slope θ_(s) going through the point x_(s),y_(s) is extrapolated to the point that is closest to the nextmeasurement point.

The system noise has a Gaussian distribution with zero mean and acovariance matrix based on independent perturbations in curvature,orientation and lateral offset (movement in a direction normal to theboundary). A Brownian model is used, where the standard deviations ofperturbations in curvature, orientation and lateral offset areproportional to the square root of the arc length of the extrapolatedcircle of the previous paragraph. The accumulated effect of curvatureperturbation on orientation and lateral offset is also modelled,resulting in the following covariance matrix:$Q = {{\sigma_{k}^{2}\begin{pmatrix}s & {\frac{1}{2}s^{2}} & {\frac{1}{6}s^{3}} \\{\frac{1}{2}s^{2}} & {\frac{1}{3}s^{3}} & {\frac{1}{8}s^{4}} \\{\frac{1}{6}s^{3}} & {\frac{1}{8}s^{4}} & {\frac{1}{20}s^{5}}\end{pmatrix}} + {s\begin{pmatrix}0 & 0 & 0 \\0 & \sigma_{s\quad\theta}^{2} & 0 \\0 & 0 & \sigma_{sy}^{2}\end{pmatrix}}}$with respect to the system curvature κ, system slope θ and lateraloffset respectively, where s is the arc length (the circular distancebetween the previous point and the estimate of the next point). Theconstants σ_(κ), σ_(sθ), σ_(sy) define the average roughness of thenuclear boundary, and depend on the type of cell being analysed.

The measurement transition matrix H maps the system parameters to themeasurement parameters in the natural way: ${\begin{pmatrix}x_{m} \\y_{m} \\\theta_{m}\end{pmatrix} = {H\begin{pmatrix}x_{s} \\y_{s} \\\theta_{s} \\\kappa\end{pmatrix}}},{H = \begin{pmatrix}1 & 0 & 0 & 0 \\0 & 1 & 0 & 0 \\0 & 0 & 1 & 0\end{pmatrix}}$

The measurement noise is based on independently Gaussian distributedperturbations of slope and lateral offset, resulting in the followingcovariance matrix with respect to measurement slope and lateral offsetrespectively: $R = \begin{pmatrix}\sigma_{m\quad\theta}^{2} & 0 \\0 & \sigma_{my}^{2}\end{pmatrix}$

The constants σ_(mθ), σ_(my) define the average smoothness of thenuclear boundary, and depend on the type of cell being analysed.

The following constants are used to define the clutter model:

-   -   ρ=Clutter density=average number of edges per unit area that are        not associated with the nuclear boundary.    -   P_(D)=Probability that the edge feature is associated with the        true nuclear boundary.

These constants depend on the clutter present in the image, both itsedge strength relative to the nuclear boundary and its average spatialdensity.

For a nucleus with average radius r the following parameters of theabove model are used (all units are image pixels):

-   -   σ_(κ)=r^(−3/2)    -   σ_(my)=3    -   σ_(mθ)=0.3 radians    -   σ_(xy)=0.9    -   σ_(xθ)=0.09 radians    -   ρ=0.01    -   P_(D)=0.8

The initial value for the system covariance matrix M is given by:${M_{0} = \begin{pmatrix}\sigma_{m\quad y}^{2} & 0 & 0 \\0 & \sigma_{m\quad\theta}^{2} & 0 \\0 & 0 & \sigma_{k0}^{2}\end{pmatrix}},{{{where}\quad\sigma_{k0}} = {1/r}}$

The PDAF estimator is now applied sequentially as follows. The matricesH_(k), Q_(k) and R_(k) are all constant in this application (as definedabove). The following expressions are those known from the Bar-Shalomand Fortmann reference quoted above.

Update

Innovation${{\underset{\_}{v^{\prime}}}_{k} = {{\sum\limits_{i = 1}^{N_{k}}{\beta_{ki}\quad{\underset{\_}{v}}_{ki}\quad{where}\quad{\underset{\_}{v}}_{ki}}} = {{\underset{\_}{z}}_{ki} - {H_{k}{\underset{\_}{\overset{\_}{x}}}_{k}}}}}\quad$Kalman Gain MatrixK _(k) =M _(k) H ^(T) _(k) S ⁻¹ _(k) where S _(k) =H _(k) M _(k) H ^(T)_(k) +R _(k)Beta Weights $\beta_{ki} = \left\{ {{{\begin{matrix}{e_{ki}/\left\lbrack {b + {\sum\limits_{j = 1}^{N_{k}}e_{kj}}} \right\rbrack} & {{{for}\quad i} \neq 0} \\{b/\left\lbrack {b + {\sum\limits_{j = 1}^{N_{k}}e_{kj}}} \right\rbrack} & {{{for}\quad i} = 0}\end{matrix}{where}\quad e_{ki}} = {{{\exp\left( {{- \frac{1}{2}}{\underset{\_}{v}}_{ki}^{T}S_{k}^{- 1}{\underset{\_}{v}}_{ki}} \right)}\quad\quad{for}\quad i} \neq 0}},{{{and}\quad b} = {{\rho\left( {1 - P_{D}} \right)}{\sqrt{{2\pi\quad S_{k}}}/P_{D}^{2}}}}} \right.$State Estimate Update{circumflex over (x)} _(k) ={overscore (x)} _(k) +K _(k) v′ _(k)Error Covariance Update$P_{k} = {{\beta_{k0}M_{k}} + {\left( {1 - \beta_{k0}} \right)P^{*}} + {{K_{k}\left\lbrack {\left( {\sum\limits_{i = 1}^{N_{k}}{\beta_{ki}\quad{\underset{\_}{v}}_{ki}\quad{\underset{\_}{v}}_{ki}^{T}}} \right) - {{\underset{\_}{v^{\prime}}}_{k}{\underset{\_}{v^{\prime}}}_{k}^{T}}} \right\rbrack}K_{k}^{T}}}$where  P^(*) = [I − K_(k)H_(k)]M_(k)PredictionState Estimate Extrapolation{overscore (x)} _(k+1)=Φ_(k) {circumflex over (x)} _(k)Error Covariance ExtrapolationM _(k+1)=Φ_(k) P _(k)Φ^(T) _(k) +Q _(k)

This process continues until the entire contour is traversed (i.e.returned to the starting point). This is then repeated around thecontour for a second time (using the final conditions from the firstpass as starting conditions for the second pass); this ensures that thefinal estimate of the smoothed contour is independent of the assumedinitial conditions.

The curvature of the smoothed contour derived from the PDAF tracker isnow used to find a measure of the degree of non-convexity of theobject's boundary. Firstly, the sign of the curvature is set so that itis positive where the boundary is locally convex and negative wherelocally concave (as viewed from outside the boundary). All positivevalues of curvature are then set to zero, leaving non-zero values onlyat locations where the boundary is locally concave. A graph of curvature(Y-axis) against perimeter arc length i.e. distance along the boundary(X-axis) is plotted, then the line integral of curvature with respect toarc length is computed. The absolute value of this integral is taken toproduce a non-negative result. The final result is a dimensionlessquantity giving an indication of overall non-convexity, called the“negative curvature area”. Objects which are almost completely convex,in this case whose negative curvature area is less than 0.2, are thenrejected.

The output from this process is a set of boundary measurements for eachobject, namely negative curvature area, and a more precise estimate ofarea.

Grey-Level Normalisation

Next a normalisation process 25 is carried out to allow for differencesin overall brightness between different slides. For each remainingobject, the mean grey level of the pixels enclosed within (but not on)the high-resolution contour found in stage 23 is calculated. Thestatistics used for normalisation are the overall mean and standarddeviation of the grey levels of the regions obtained from stage 22. Eachobject's grey level is then normalised (by subtracting this mean anddividing by this standard deviation). The output is a statistic for eachassumed nucleus.

Second Clutter Rejection

The next process 26 involves a second stage of classification andclutter rejection based on the Fisher classifier to discriminate betweenobjects representing mitotic and non-mitotic nuclei. The Fisherclassifier is a known statistical classification method described forexample in Section 4.3 of “Statistical Pattel Recognition” by Andrew R.Webb, Arnold Press, 1999, and is preferred for this stage of the processdue to its robustness against overtraining.

In this case the Fisher classifier uses a set of information about eachobject that has been derived by analysis as described above. Eachinformation set is a feature vector, that is an ordered list of realnumbers each describing some aspect of the object; each component numberis denoted an object feature. The purpose of the classificationalgorithm is to discriminate between two classes of object based on theinformation contained in their feature vectors. The output of thealgorithm is a set of numbers, one for each object, indicating thelikelihood that the nucleus which it represents is a member of one ofthe two chosen classes (in this case the classes are mitotic andnon-mitotic).

For a given feature vector x, the standard implementation of the Fisherclassifier output is defined as:$F = {\sum\limits_{k = 1}^{n}{a_{k}x_{k}}}$where x=[x₁, . . . x_(k)] is the feature vector. In this embodiment,this definition has been extended to use non-linear functions of thefeature vector, namely: $F = {\sum\limits_{k = 1}^{n}{a_{k}{g_{k}(x)}}}$where a_(k) are prescribed real numbers and g_(k) are prescribedfunctions of the feature vector x. These functions and variables arechosen to give the lowest number of misclassifications for objects withknown class.

The components of the feature vector x are mean grey level (computed instage 25) and negative curvature area (computed in stage 24). Inprinciple area could also be used, since smaller mitotic cell nucleitend to be darker and less concave than larger ones. In this case aquadratic set of basis functions (g_(k)) are used, so that the Fisherclassifier value is given by:F=a ₁ G+a ₂ G ² +a ₃ C+ a ₄ GC+ a ₅ C ² +a ₆where G is the normalised grey-level, C is the negative curvature area,and the coefficients a₁ are derived from the training stage referred tobelow.

The coefficients and decision boundary for the Fisher classifier areobtained by training the classifier on a large number of example slidesprovided by a histopathologist where accurate ground truth (sets ofmitotic and non-mitotic cells) is also available. The training stageresults in a classifier boundary which minimises the total number ofmisclassifications, i.e. both false negatives (missed mitotic cells) andfalse positives (falsely-detected non-mitotic cells). In the preferredembodiment the resulting coefficients ai that have been derived fromthis training stage are [−0.87431, 0.10205, 0.84614, −0.18744, −0.04954,−5.56334]. FIG. 5 illustrates the classifier together with the data onwhich it was trained, where plusses indicate mitotic cells and crossesindicate non-mitotic cells.

Mitosis Count

Stage 27 counts the number of objects deemed to represent the nuclei ofmitotic cells, that is to say only those objects whose values exceed agiven threshold in the output of the Fisher classifier. The preferredcriterion is F>0 (illustrated as the decision boundary in FIG. 5), setto give the optimum trade-off between missed mitotic cells andfalsely-detected non-mitotic cells. The number of objects whoseclassifier value exceeds this threshold defines the mitotic count forthat tile. The count for the ten tiles analysed is aggregated, and canbe converted into an indication of cancer grade in accordance with atable as previously described in connection with the existing manualprocedure.

It will be appreciated that the coding of a computer program toimplement all of the processing stages described above for the preferredembodiment of the invention can be achieved by a skilled programmer inaccordance with conventional techniques. Such a program and code willtherefore not be described further.

1. A method for the automated analysis of a digital image comprising anarray of pixels, including the steps of: (a) identifying the locationsof objects within the image which have specified intensity and sizecharacteristics; (b) defining regions of specified extent within theimage which contain respective said objects; (c) deriving from the datawithin respective said regions one or more respective closed contourscomprising points of equal intensities; and (d) estimating the curvatureof at least one respective said contour within respective said regionsat least to produce a measure of any concavity thereof.
 2. A methodaccording to claim 1 wherein step (a) comprises the application of aradially-symmetric difference filter with zero mean.
 3. A methodaccording to claim 2 wherein the image is filtered at a plurality ofresolutions of increasing scale.
 4. A method according to claim 2wherein said locations are identified in accordance with the locationsof respective local extrema in the output of said filter.
 5. A methodaccording to claim 4 including the step of sorting, in order ofintensity, local extrema in the output of said filter and selecting forfurther analysis only those objects which correspond to a specifiedproportion of said extrema in such order.
 6. A method according to claim1 further comprising, following step (a): selecting an intensitythreshold related to the mean intensity of pixels within the image inneighbourhoods of said locations; creating a binary image according towhether pixels in the first-mentioned image are above or below saidthreshold; identifying regions in the binary image composed of connectedpixels which are below said threshold in the first-mentioned image; andrejecting from further analysis those objects which correspond to suchregions in the binary image which fall below a specified size orthickness.
 7. A method according to claim 1 wherein step (c) comprises,for respective said regions, deriving respective first and second saidcontours having respectively lower and higher resolutions, determiningwhether the sizes and locations of said first and second contours areconsistent within specified criteria and, if so consistent, selectingsaid second contour for step (d).
 8. A method according to claim 7wherein, for respective said regions, the first said contour is derivedby: seeking within the region one or more contours of respectivespecified intensities; determining whether the or each such contour is aclosed contour and meets specified location, size and/or intensityorientation criteria; and if more than one such contour is a closedcontour and meets such criteria, selecting from the same the contour ofthe lowest intensity.
 9. A method according to claim 8 wherein saidspecified intensities are no greater than that which corresponds to thecontour of highest edge strength within the respective region.
 10. Amethod according to claim 9 wherein step (a) comprises the applicationof a radially-svmmetric difference filter with zero mean and said firstcontour is derived by seeking one or more contours in the output of saidfilter for the respective region and said specified intensities are nogreater than the zero level in such output.
 11. A method according toclaim 8 wherein, for respective said regions, the second said contour isderived by: seeking within the region a plurality of contours ofrespective specified intensities ranging between the lowest and highestintensities within the region; determining whether each such contour isa closed contour and meets specified location, size and/or intensityorientation criteria; and if more than one such contour is a closedcontour and meets such criteria, selecting from the same the contourhaving the highest edge strength.
 12. A method according to claim 1wherein step (d) includes the application of a Probability DensityAssociation Filter to respective said contours.
 13. A method accordingto claim 1 wherein step (d) comprises, for respective said contours:measuring the curvature of the contour at a plurality of points aroundthe contour, convexity and concavity being of opposite sign; settingconvex values of such curvature to zero; plotting resultant values ofcurvature at said points against a measure of the distance of therespective point along the contour; and computing as said measure ofconcavity the line integral of such plot.
 14. A method according toclaim 1 further comprising the step of: (e) classifying objects into oneof at least two classes in accordance with a function of said measure ofconcavity of a contour corresponding to the respective object and ameasure of the mean intensity of the respective object.
 15. A methodaccording to claim 14 wherein step (e) is performed by use of a Fisherclassifier.
 16. A method according to claim 14 wherein the intensitiesof respective objects are normalised prior to step (e).
 17. A methodaccording to claim 14 further comprising the step of: (f) counting thenumber of objects classified into a specified one of said classes.
 18. Amethod according to claim 1 wherein the image is of a histological orcytology specimen or of a soil sample.
 19. (canceled)
 20. A methodaccording to claim 17, wherein the image is of a section of breasttissue and said specified class is identified as the class of mitoticepithelial cell nuclei.
 21. (canceled)
 22. A method for the automatedidentification of mitotic activity from a digital image of ahistological specimen, including the steps of: (a) identifying thelocations of objects within the image which have specified intensity andsize characteristics associated with epithelial cell nuclei; (b)defining regions of specified extent within the image which containrespective said objects; (c) deriving from the data within respectivesaid regions one or more respective closed contours comprising points ofequal intensities; (d) estimating the curvature of at least onerespective said contour within respective said regions at least toproduce a measure of any concavity thereof; and (e) classifying objectsas representing mitotic cell nuclei as a function of at least saidmeasure of concavity of a contour corresponding to the respectiveobject. 23-24. (canceled)
 25. A computer program comprising instructionsto cause a computer to execute a method according to claim 1.